#ifndef _SP_SOLVER_K_H_
#define _SP_SOLVER_K_H_
#include"help.h"
#include"orbital.h"
#include"hfsolver_ho.h"
#include"system_hf.h"
#include"constants.h"
#include<limits>
#include<vector>
#include<eigen3/Eigen/Eigen>//matrix library
#include<complex>

using std::abs;
using std::numeric_limits;
using std::stable_sort;
using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

using Eigen::ComplexEigenSolver;

using std::vector;
using std::complex;

typedef complex<double> complexd;

//typedef double KDataType;
typedef complexd KDataType;
typedef Matrix<KDataType,Dynamic,Dynamic> KMat;
typedef Matrix<KDataType,Dynamic,1> KVec;
//typedef SelfAdjointEigenSolver<KMat> EigenSolver;
typedef ComplexEigenSolver<KMat> EigenSolver;


/// SP. Orbital wavefunction in K representation
class SP_K:public Orbital_ljjtz
{
 public:
  SP_K(){}
 SP_K(int _l,int _jj,int _tz,int _channel,KDataType _e,KVec _wavefun):Orbital_ljjtz(_l,_jj,_tz),channel(_channel),e(_e),wavefun(_wavefun){}
  int channel;
  KDataType e;
  KVec wavefun;
  int type;
  bool operator < (const SP_K & sp) const
  {
    //in case that energy is too close
    if( abs(e-sp.e) > numeric_limits<double>::epsilon() )
      {
    	return (e < sp.e);
      }
    else
      {
    	return false;
      }
  }

  //for basis transform, overlap between this and HO. orbitals   
  int GroupIndex;
  vector<KDataType> Coeff;
};

class SPChannel_K
{
 public:
  SPChannel_K(){}
  SPChannel_K(int _l,int _jj,int _tz,int _numPoints):l(_l),jj(_jj),tz(_tz),numPoints(_numPoints),k(_numPoints),wk(_numPoints){}
  int l,jj,tz;

  int numPoints;///<number of quadration points
  ///quadration points and weights
  vector<KDataType> k;
  vector<KDataType> wk;
};

class SPSolver_K
{
 public:
 SPSolver_K():hfsolver(SystemHF.hfsolver),rdim(80),rx(rdim),rw(rdim){}
 SPSolver_K(System_Tz_HO * Ptr,int _A,int _Z):rdim(80),rx(rdim),rw(rdim),pSystem_Tz(Ptr),A(_A),Z(_Z),realA(_A),SystemHF(Ptr,_A,_Z),hfsolver(SystemHF.hfsolver)
    {
      gauleg(0,10,rx,rw);
    }

  void setup(System_Tz_HO * Ptr,int _A,int _Z)
  {
    pSystem_Tz=Ptr;    
    A=_A;
    Z=_Z;
    realA=_A;
    SystemHF.setup(Ptr,_A,_Z);
  }
  void setAZ(int _A,int _Z,int _realA=-1)
  {
    if(_realA==-1) realA=_A;
    else realA=_realA;
    A=_A;
    Z=_Z;
    SystemHF.setAZ(A,Z,realA);
  }
  
  void setup();

  KDataType getVkk(int channel,int a,int b);///< return potential V_{ljtz}(k[a],k[b]),ljtz is quantum numbers of this channel

  ///get HF. orbital in k space at point k of the i'th orbital
  KDataType hf_k(int i,KDataType k);
  ///get sum_k {<k1 h|k1 \dot k2|k2 h>}
  KDataType getHFk1k2(KDataType k1,KDataType k2,int group);

  void setVkkMat(int channel,KMat &h);
  void setVkkMatfromHF(int channel,KMat &h);

  double getVrrfromHF(double r1,double r2,int group);
  double getVrrfromHF(int channel,double r1,double r2);
  KDataType getVkkfromHF(int channel,int a,int b);
  KDataType getVkkfromVrr(int channel,int a,int b);
  KDataType getVrrfromVkk(int channel,double r1,double r2);
  void printVrr(int channel);
  void printHFVrr(int numpoints,double step=0.1);//print HF potential


  
  const int rdim;//quadrature for bessel-fourier transform
  vector<double> rx,rw;
  double getVrr(int channel,double r1,double r2);
  
  
  void solve();

  void setupOverlap(int num);///< set up overlap with HO. orbitals, the first "num" orbitals is set up

  ///print total 'num' orbitals, if 'channel' is offered, only print orbitals of that channel. default print all orbitals of all channels
  void printOrbitals(int num=-1,int channel=-1);
  void printVkk(int channel);



  KDataType get1B(int bra,int ket) const;
  KDataType get1B_p2(int bra,int ket) const;
  
  vector<SPChannel_K> Channels;
  
  vector<SP_K> Orbitals;


  void setType();///< decide the orbital type, only for rotation contour.


  
  
  int A,Z;
  int realA;
  int FermiSurface;
  KDataType Energy;
  
  System_Tz_HO * pSystem_Tz;
  const SHFSolver_HO & hfsolver;

  System_Tz_HF SystemHF;
};

#endif
